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Abstract 

Using  conventional  particle  tracking  techniques  to  predict  the  dynamics  of  spray  flows  can  be  prohibitively 
expensive,  requiring  large  computation  times  and  significant  data  storage.  Moreover,  because  of  the  discontinuous 
nature  of  the  spray  drops,  data  from  a  simulation  of  the  flow  does  not  produce  smooth  statistics  unless  the  results 
from  many  simulations  have  been  averaged.  Recently,  a  new  model  was  developed  that  computes  spray  statistics 
directly,  without  simulating  the  flow.  In  this  paper,  the  model  is  extended  to  include  the  effects  of  nonlinear 
momentum  exchange  between  the  phases.  The  approach  was  tested  on  a  quasi-one-dimensional  flow  geometry. 
The  results  are  compared  with  a  Lagrangian  simulation  and  demonstrate  good  agreement. 


Introduction 

Spray  flows  play  an  important  role  in  many 
industrial  processes  today.  It  is  therefore  vital  that  we 
have  a  full  understanding  of  the  physics  of  these  flows 
and  are  able  to  accurately  predict  both  the  distribution 
and  behavior  of  the  drops  and  the  dynamics  of  the  gas 
phase.  Failure  to  adequately  model  these  flows  could 
lead  to  inefficient  fuel  mixing,  non-uniform  industrial 
coatings,  or  poorly  distributed  agricultural  sprays. 
Additionally,  since  one  may  wish  to  make  small 
changes  to  an  engine  or  injector  design  to  determine  the 
effects  on  the  two-phase  mixture,  it  is  of  significant 
interest  to  reduce  the  time  required  to  predict  these 
spray  flows. 

Perhaps  the  most  common  approach  to  studying 
spray  flows  has  been  the  Lagrangian-Eulerian,  or 
particle-tracking  method.1,2  3,4  Drops  are  stochastically 
injected  into  the  gas  phase  and  their  trajectories  are 
determined  by  integrating  a  Lagrangian  equation  of 
motion.  The  gas  phase  is  typically  modeled  with  the 
time-dependent  RANS  equations  with  a  suitable 
turbulence  model  and  exchange  terms. 

While  particle-tracking  methods  have  provided 
useful  information  in  many  applications,  they  have 
some  potentially  significant  drawbacks.  For  instance, 
data  arising  from  a  simulation  must  be  post-processed. 
If  the  quantity  of  interest  is  the  mean  number  density  in 
each  grid  cell,  this  may  not  pose  a  problem.  However, 
if  we  are  interested  in  more  detailed  statistics,  the 


required  computational  time  increases  because  a 
sufficient  number  of  drops  must  pass  through  the  cell  to 
provide  a  data  set  large  enough  for  a  meaningful 
average.  As  the  quantity  of  interest  becomes  more 
specific,  the  necessary  computation  time  becomes  more 
prohibitive. 

An  alternative  approach  which  does  not  involve 
simulation  is  to  compute  the  evolution  of  a  probability 
density  function  (PDF)  describing  the  drops.  Williams5 
was  the  first  to  derive  a  transport  equation  for  a  droplet 
PDF,  called  the  spray  equation,  analogous  to 
Boltzmann’s  equation  for  molecules.  Direct  solutions 
of  the  spray  equation  have  been  attempted,  ’  ’  ’  but 
only  with  limited  success.  The  high-dimensionality  of 
the  equation  and  current  computational  resources  limit 
the  numerical  resolution  attainable  within  the  phase 
space. 

Recently,  a  Maximum  Entropy  Moment  Closure 
(MEMC)  model  was  described  that  gives  a  complete 
description  of  a  spray  flow  by  computing  the  evolution 
of  its  PDF  along  with  the  gas  flow  in  which  it  is 
embedded.11  It  was  shown  that  both  general  and 
detailed  statistics  about  the  spray  drops  and  the  solution 
of  the  gas  phase  could  be  attained  without  the  need  to 
average  simulation  data  and  in  significantly  less  time 
than  might  be  required  in  a  typical  particle-tracking 
method.  The  purpose  of  this  paper  is  to  incorporate  a 
nonlinear  drag  law  into  the  model  and  compare  test 
results  to  a  Lagrangian  simulation. 
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Governing  Equations 

The  equations  that  describe  the  evolution  of  our 
system  make  use  of  ensemble  averaging  which  is 
independent  of  any  time  or  length  scales.  They  also 
restrict  us  to  single-point  statistics.12  Each  drop  is 
described  by  a  relatively  small  set  of  characteristics 
which  contains  the  information  of  interest.  This  set  of 
characteristics  describes  a  state  vector  a  whose 
elements  define  the  axes  of  a  hyperspace,  called  a- 
space,  through  which  the  drop  can  move. 

Because  we  are  dealing  with  point  particles, 
volume  displacement  effects  are  negligible,  and  we  are 
able  to  treat  any  interphase  exchange  as  if  the  drop  were 
alone  in  a  locally  uniform  gas  field. 

With  these  assumptions,  a  probabilistic  description 
of  the  spray  can  be  developed,12  defining  X(x;t)f  as 
the  probability  density  of  finding  a  single  drop  in  the 
vicinity  of  a  point  (x,/1).  It  is  also  the  expected 
number  density  of  drops  at  x  and  t.  This  dual 
interpretation  is  important  since  we  are  generally  more 
interested  in  the  expected  number  density  rather  than 
the  probability  density.  The  probabilistic  description 
also  defines  the  function  f(a'\xj),  the  probability 
density  of  finding  a  particular  diameter  and  velocity. 
This  is  conditioned  on  a  drop  existing  at  x  and  /,  and 
a  -space  is  the  subset  of  a  -space  excluding  the  spatial 
coordinate.  By  combining  these  two  functions,  we 
define  the  unnormalized  single-particle  probability 
density  function  for  the  spray 

F(a\t)da  -  X(x  ;t)f(a;x9t)dadV  (1) 

F(a;t)  is  unnormalized  because  X[x9t)9  also  an 
unnormalized  PDF,  determines  how  much  spray  is 
present  at  x  and  t. 

If  a  =  (x9<p,v9T) ,  a  transport  equation  for  F(a;t) 
can  now  be  derived.13’14,15 

^  +  Vx.(Fv)+Vv.(f2 

a  a  (2) 

+^(™)+^F(/r0)=^1+*0 

Equation  (2)  is  the  spray  equation.  It  describes  the 
evolution  of  the  probability  density  function 
F(x9ty9\9T9t)  through  joint  physical,  diameter, 


f  For  clarity  when  denoting  a  PDF,  those  arguments  over  which  the 
function  is  a  density  will  be  listed  first,  followed  by  a  semi-colon, 
followed  by  arguments  that  are  parameters  of  the  function. 


velocity,  and  temperature  space  and  includes  source 
terms  accounting  for  binary  collision  S2,  unary  breakup 
Sy ,  and  zero-body  events  SQ  such  as  nucleation  and 
complete  vaporization.  The  spray  equation  has  not 
generally  been  viewed  as  a  practical  way  of  predicting 
spray  flows.  It  is  an  unusual  evolution  equation  in  the 
sense  that  the  quantity  of  interest  is  not  only  being 
transported  through  physical  space,  but  also  through 
diameter,  velocity,  and  temperature  space.  A  numerical 
solution  requires  a  full  discretization  of  this  hyperspace, 
making  a  direct  solution  difficult. 

Considering  a  non-vaporizing,  isothermal  spray 
without  collisions,  breakup,  or  zero-body  events,  the 
spray  equation  reduces  to 

— +V  •(Fv)  +  Vv.(F5)  =  0  (3) 

dt 


The  diameter  axis  is  then  discretized  into  equal-sized 
bins,  taking  the  approach  of  Tambour.16  By  splitting 
and  integrating  equation  (3),11  a  set  of  velocity  moment 
transport  equations  within  each  diameter  bin  can  be 
derived  as 


+  V(A„  (vv,)„  )  =  Zj-W.  W 


3(A»(V<)„) 


dt 

3(A»(V>)J 

dt 


(5) 


dt 


+v. <7> 


3U,(v,vr)  )  ,  ,  ,  x 

V— +M^v*v-A) 

=^((a^l+(a‘vA) 


(8) 


To  determine  the  droplet  diameter  distribution,  an 
evolution  equation  for  the  binned  droplet  number 
density  is  solved 


(9) 


(16) 


Xn  is  the  probability  of  finding  a  drop  in  the  vicinity  of 
some  spatial  point  at  some  instant  in  time  within  the  rih 
diameter  bin.  The  overall  droplet  number  density  is 
given  by 


kv*)„  =-~((v*}„  -(«*v*>„ 

1  n 

Quasi-One-Dimensional  Spray  Problem 


a  =  £a„  (10) 

n 


The  terms  on  the  right  sides  of  equations  (4) 
through  (8)  account  for  the  momentum  exchange 
between  the  gas  and  liquid  phases.  The  acceleration 
term  a  is  given  by 


a  = 


where 


giving 


CD  =— (l+0.15Re0687  ) 
0  Re '  ' 


(ID 


(12) 


18p,v, 


p,f 


r\\-u\<p  ^ 

0.687  “ 

1  +  0.15 

V, 

_ 

k  ) 

_ 

(v— w)  (13) 


Defining  g-v-u  as  the  slip  velocity,  the  x- 
component  of  the  expected  acceleration  in  the  n,h 
diameter  bin  is 


-(r3i3-r3,3)(ifr^,>„ 


(14) 


where 


T  = 


_  r  IT  geo.  n 


,8PA- 


(15) 


is  the  representative  droplet  relaxation  time  of  the  n'" 
diameter  bin  and  0,  and  02  are  the  starting  and  ending 
diameters  of  the  bin,  and  n  =  ^0,02  .  Similarly, 


Because  we  are  interested  in  a  spray  problem  that 
is  understandable,  physically  interesting,  and  tractable 
while  at  the  same  time  being  able  to  evaluate  the  spray 
model,  we  make  use  of  a  quasi-one-dimensional  flow. 
Figure  1  shows  the  geometry  for  this  problem.  At 
x  =  0  there  is  a  two-dimensional  array  of  spray 
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Figure  1  Geometry  for  quasi- 1 D  spray  flow 

injectors  that  extends  infinitely  in  the  y-direction.  Each 
injector  delivers  a  distribution  of  drops  into  an 
incompressible  gas,  which  initially  has  a  low  level  of 
turbulence.  The  array  and  each  injector  on  it  extend 
infinitely  in  the  z-direction,  resulting  in  a  series  of 
wedge  sprays  of  spherical  drops. 

Now  imagine  an  ensemble  of  these  arrays/ 
Because  we  wish  to  develop  a  quasi-one-dimensional 
flow,  we  restrict  the  PDF  for  the  drop  velocities  at  each 
injector  to  have  a  mean  transverse  (y-direction)  velocity 
of  zero.  If  there  is  an  infinite  number  of  realizations  in 
the  ensemble,  then  there  is  no  preferential  location 
along  the  transverse  axis  since  there  is  equal  probability 
for  a  drop  to  have  a  positive  y-velocity  of  some 
magnitude  as  there  is  for  a  drop  to  have  a  negative  y- 
velocity  of  the  same  magnitude.  Similarly,  we  have  a 
zero-mean  gas  velocity  component  in  the  y-direction, 
though  fluctuations  do  occur.  Because  there  is  no 
preferential  location  along  the  transverse  axis,  there  are 
no  variations  along  that  axis  across  the  ensemble,  and 
therefore,  derivatives  of  averaged  quantities  in  the  y- 


+ 

+  To  ensure  a  quasi-one-dimensional  geometry,  two  approaches  can 
be  taken  with  the  injectors.  First,  in  each  realization  of  the  ensemble, 
the  injectors  can  be  placed  on  the  array  at  random.  Alternatively,  the 
injectors  could  be  equally  spaced  on  each  array,  but  shifted  slightly 
from  one  realization  to  another.  In  either  case,  the  flow  will  be  quasi- 
one-dimensional  across  the  ensemble. 


direction  are  zero.  So,  the  three  primary  constraints  for 
this  problem  are 

(»,)  =  0,  (Vj.)  =  0,  m.'  (17) 

Also,  so  that  no  mean  flow  develops  in  the  transverse 
direction,  there  can  be  no  correlation  between  the  x -  and 
^-components  of  velocity  in  either  phase.  Thus,  all 
cross-component  and  cross-component/cross-phase 
moments  are  zero.  Substituting  these  constraints  into 
the  transport  equations  yields 


3A, 

dt 


f-+ Vx  -(A„  (vx }/7)~0 


(18) 


mw.)  AW).) 


dt 

A 

T. 

8. 6  A. 


dx 


“<«.))  09) 


Aft. 


dt 


dx 


=-7l((vv)„-Kv^U 


(20) 


17.2A„ 


p,v"n(r 


-o(i*r*.v.)a 


s(^(v?)j,a(^(v^).) 

3/  3x 

=  ~~L({s/l)„~{uyvy)„)  (21) 

Ln 

pX3  13  (r313  -  r 3,3  )(|gr7  s,v...  \ 


If  constraints  (17)  are  applied  to  the  incompressible 
RANS  equations,  we  obtain 


0 


(22) 


dt  dx  dx  (23) 

+3itvxX(tp  (vr -«r)(l  +  O.I5  Re06*7 )) 


Equation  (22)  says  that  the  mean  gas  velocity  is  a 
constant  in  space.  If  we  specify  the  inlet  mean  gas 
velocity  as  a  constant  in  time,  then  the  mean  gas 
velocity  throughout  the  domain  is  constant  in  both 
space  and  time.  This  makes  solving  the  gas  momentum 
equation  unnecessary  unless  one  is  interested  in  the 
pressure. 

Closure 

Several  models  are  used  to  close  the  above  set  of 
equations.  First,  a  maximum  entropy  approach15,  l6,  17 
is  taken  to  close  the  higher-order  moments,  such  as  the 
third-order  velocity  moment  in  the  advection  term  in 
equation  (20).  The  maximum  entropy  principle  states 
that  of  all  possible  PDFs  consistent  with  a  set  of  given 
(moment)  constraints,  the  one  PDF  that  is  most 
unbiased  is  the  one  with  the  most  entropy.  All  known 
information  is  included  in  the  form  of  the  constraints, 
with  no  extraneous  information  that  might  tend  to  bias 
the  PDF.  We  can  then  integrate  the  maximum  entropy 
PDF  to  obtain  values  for  the  higher-order  moments. 

Additional  modeling  is  required  for  the  cross¬ 
phase/same-component  moments  that  appear  in  the 
equations.  In  the  derivation  of  the  spray  equation,  the 
gas  velocity  was  not  included  as  part  of  the  set  of 
characteristics  that  describes  a  drop.  This  amounts  to 
an  independence  assumption  between  the  two  phases 
which  is  non-physical.  We  must  therefore  develop  a 
model  for  these  moments.  This  will  be  done  by 
choosing  a  functional  form  for  the  correlation 
coefficient.  The  correlation  coefficient  is  defined  as 


Corrx  = 


(wrv, 


)*  (wx)(vx)„ 


(24) 


where  ou  and  crv  are  the  gas  and  drop  velocity 

variances,  respectively,  in  the  nh  diameter  bin.  In 
general,  the  correlation  coefficient  can  have  any  value 
between  negative  one  and  positive  one,  inclusive.  The 
physics  of  this  problem,  however,  dictate  that  the 
coefficient  must  be  bounded  between  zero  and  positive 
one.  Physically,  this  is  because  a  positive  increase  in 
the  drop  velocity  will  produce  a  positive  increase  in  the 
gas  velocity. 

Perfect  correlation  between  the  velocities  of  the 
two  phases  implies  that  their  means  and  variances  are 
identical.  Any  differences  would  reduce  the  degree  of 
correlation,  and  since  those  differences  can  be  large,  we 
assume  that  there  is  zero  correlation  if  either  the 
difference  of  the  means  or  the  differences  of  the 
variances  are  infinite.  We  also  know  that  the 
correlation  depends  upon  the  drop  size.  The  velocities 
of  small  drops  ten  to  be  more  correlated  with  the  gas 
velocity  because  of  their  small  relaxation  times. 


(29) 


whereas  the  velocities  of  larger  drops,  which  take 
longer  to  respond  to  fluctuations  in  the  gas  phase,  tend 
to  be  less  correlated.  All  of  these  considerations  are 
accounted  for  if  the  correlation  coefficient  is  modeled 
as  an  exponential  function 


/ 

Corrx  =  exp 

v 


\ 


) 


(25) 


where  the  ratio  of  time  scales  represents  a  droplet 
Stokes’  number.  Ca  is  an  adjustable  constant  that  must 
be  determined  from  experiment,  xn  is  the  drop 
relaxation  time,  and  xg  is  the  gas-phase  time  scale. 
Short  drop  relaxation  times  will  have  velocities  that  are 
highly  correlated  with  the  gas  velocity.  Thus,  drops  in 
the  limit  of  having  no  mass  will  be  perfectly  correlated. 
Large  drops  with  long  relaxation  times  will  not  be  as 
well  .correlated.  The  gas  time  scale  is  taken  to  be  the 
shorter  of  the  mean  residence  time  of  a  drop  in  a 
turbulent  eddy  with  size  equal  to  the  turbulent  length 
scale,  or  the  turbulent  time  scale.  If  the  residence  time 
within  an  eddy  is  short,  that  means  that  there  is  strong 
slip  between  the  drops  and  the  gas,  resulting  in  a  low 
velocity  correlation.  On  the  other  hand,  if  the  mean 
residence  time  is  long,  then  the  correlation  is  based 
upon  the  eddy  lifetime,  as  defined  by  the  turbulent  time 
scale. 

To  find  the  turbulent  time  and  length  scales,  and  to 
close  the  Reynolds  stress  term  in  equation  (23),  we 
introduce  a  k—e  model  that  has  been  modified  to 
account  for  the  presence  of  the  drops.  Applying 
constraints  (17),  the  governing  equations  for  the  model 
are 


dt  a* 


V..  +-T- 


*  4 


(26) 


de  |  3(g(«)) 

dt  dx 

_-Ctf2g  +  Cf3y,  ,  3 


T  dx 

where  T  is  the  turbulence  time  scale. 


v,  )de 
v  +— - 

*  U* 


(27) 


v,  = 


C„*2 


(28) 


and 


C„  =0.09  Ce2=  1.92 
C„=  1.5  4=1-0  4=1-3 

where  the  constant  values  are  the  same  as  in  Amsden  et. 
al?  Notice  that  the  terms  containing  the  production  of 
turbulent  kinetic  energy  are  identically  zero  under  the 
imposed  geometric  constraints.  The  term 

-^=(FX)  +  (F>;)  (30) 

is  the  negative  of  the  rate  at  which  the  turbulent  eddies 
are  doing  work  dispersing  the  drops.  If  Ws  <  0 ,  then 
the  turbulent  kinetic  energy  is  being  transferred  from 
the  gas  to  the  drops,  and  if  Ws  >  0  ,  the  drops  are  losing 
their  energy  and  transferring  it  to  the  gas  in  the  form  of 
turbulent  kinetic  energy.  F'  and  F'  are  the 

fluctuations  from  the  mean  force  in  the  x -  and  y - 
directions  acting  on  the  gas  by  the  drops.  Equation  (30) 
is  computed  by 

Ws  =-f  F(a;t)Fuda  (31) 

Jd' 

For  compactness,  we  do  not  show  the  final  expression 
here  as  it  is  rather  lengthy.  However,  we  should  note 
that  it  contains  numerous  nonlinear  terms  that  would  be 
too  expensive  to  compute  directly.  To  remedy  this, 
each  term  was  expanded  in  a  Taylor  series  thus 
reducing  the  time  required  to  compute  these  terms. 

Solution 

The  moment  equations  are  solved  using  an 
implicit,  first-order,  upwind  scheme  for  the  spatial 
derivatives  and  a  second-order  backwards  difference  for 
the  time  derivatives.  A  first-order  scheme  was  chosen 
for  the  spatial  derivatives  due  to  the  difficulties 
associated  with  solving  moment  equations.  Moment 
variables  have  relations  among  them  that  must  be 
preserved,  constraining  their  values.  For  example,  the 
relation  between  the  first  and  second  moments  of  a  PDF 
is 

(x2)-(x)2>0  (32) 

where  X  is  any  random  variable.  If  this  relation  is  not 
maintained,  then  the  numerical  solution  may  either 
become  unstable  or  give  non-physical  answers. 
Second-order  numerical  solutions  of  hyperbolic 
equations  can  cause  oscillations  in  the  solution  in 
regions  of  strong  gradients,  even  when  artificial 
dissipation  is  added.  Being  non-physical,  those 
oscillations  could  cause  the  numerical  solution  of  the 
variance  to  become  negative,  corrupting  the  remainder 


of  the  solution.  Therefore,  the  first-order  numerical 
scheme  is  presently  used. 

Results 

To  provide  a  comparison  between  the  MEMC 
model  and  more  conventional  methods,  a  Lagrangian 
simulation  of  the  quasi-one-dimensional  spray  is 
performed.  Parcels,  each  representing  100  drops,  are 
tracked  from  the  injector  along  their  trajectories. 
Because  the  mean  gas  velocity  is  a  constant  in  space 
and  time,  the  gas  momentum  and  continuity  equations 
are  not  solved,  though  the  k-e  model  is  used  to 
simulate  the  turbulence,  giving  each  drop  a  random 
“kick”  over  a  duration  equal  to  the  eddy  turnover  time 
or  the  drop  residence  time,  whichever  is  shorter.  This  is 
similar  to  what  is  done  in  other  Lagrangian  simulations 
(e.g.,  KIVA3).  The  trajectory  equations  are  solved  by  a 
fourth-order  Runge-Kutta  method.  Results  for  four 
simulations  are  averaged  and  only  shown  for  0.03  sec 
of  flow  time.  This  is  due  to  the  extensive  time  required 
to  compute  Lagrangian  statistics  at  longer  flow  times. 

Table  1  lists  the  moments  used  to  specify  the 
droplet  PDF  at  the  injector  and  the  gas  phase 


Quantity 

Value 

3.0  m/s 

crv 

3.0  m/s 

{",) 

6.0  m/s 

M 

0.0  m/s 

au 

0.3  m/s 

0.3  m/s 

Quantity 

Value 

X 

1000  drops 
cm 

<*> 

100  pirn 

<*s> 

150  Urn2 

(?) 

300  fim3 

<v,> 

30  m/s 

<v,> 

0.0  m/s 

Table  1  Injector  Conditions 

characteristics.  At  the  injector,  there  is  a  24  m/s  slip 
velocity  between  the  phases,  well  into  the  non-linear 
drag  regime.  The  diameter  PDF  is  shown  in  Figure  2. 
This  PDF  was  produced  using  the  Maximum  Entropy 
Formalism  with  the  constraints  listed  in  Table  1.  The 
velocity  PDFs  are  Gaussian,  which  is  also  the 
maximum  entropy  solution  for  the  given  moment 
constraints  in  Table  1.  The  inlet  condition  for  the 
turbulent  kinetic  energy  corresponds  to  a  RMS  value  of 
five  percent  of  the  mean  gas  velocity  and  the  rate  of 
dissipation  of  turbulent  kinetic  energy  is 
e  -  kl  100  sec  . 

Figure  3  shows  profiles  of  the  mean  number 
density  at  various  instants  in  time.  The  spray 


Figure  2  Injector  diameter  PDF 


Figure  3  Expected  drop  number  density 


propagates  through  the  domain  as  a  concentration  wave 
as  indicated  by  the  advancement  of  the  leading  edge. 
Starting  with  1000  drops/cm3  at  the  injector,  the 
expected  number  density  increases  downstream  as  the 
drops  decelerate  in  the  slower  moving  gas  until  we 
reach  the  leading  edge  where  it  drops  off  sharply. 
There  is  excellent  agreement  with  the  results  from  the 
Lagrangian  simulation  at  0.03  sec.  Notice,  however, 
the  amount  of  noise  present  in  the  simulation  data.  This 
is  indicative  of  this  type  of  calculation  where  it  is 
necessary  to  compute  discrete  data,  whereas  in  the 
MEMC  model,  we  get  a  smooth  curve  representing  the 
averaged  data. 

At  the  leading  edge,  the  mean  number  density  is 
quite  low.  When  drops  first  arrive  at  some  spatial 
location,  only  a  few  are  present,  those  that  were  large 
enough  not  to  have  been  greatly  influenced  by  drag. 
This  is  shown  in  Figure  4  where  the  expected  drop 
diameter  is  plotted  At  the  leading  edge,  there  is  a 
greater  mean  drop  diameter.  This  is  caused  by  the 
strong  slip  velocity  present  in  the  flow,  decelerating  all 
but  the  most  massive  drops.  As  we  move  upstream,  we 


MEMC 


Figure  4  Expected  drop  diameter 


encounter  the  smaller  drops  that  were  quickly 
decelerated. 

Figure  5  shows  how  the  standard  deviation  of  the 
diameter  PDF  evolves.  At  the  leading  edge,  the 
standard  deviation  is  large  due  to  the  fact  that  there  are 
so  few  drops  out  on  the  wave  front.  Those  drops 
present  are  large  due  to  the  size- velocity  correlation  that 


Figure  5  Drop  diameter  standard  deviation 


has  developed.  Figure  2  shows  that  there  is  a  wide 
range  of  large  diameters  that  those  drops  have.  It  can 
be  argued  that  the  dip  in  the  standard  deviation  curve 
behind  the  leading  edge  is  a  phenomenon  that  arises 
from  the  spatial  spreading  of  the  drops  that  were 
initially  injected,  however,  due  to  the  diffusion  in  the 
numerical  solution  at  the  leading  edge,  this  effect  may 
appear  more  significant  than  it  really  is.  Finally,  notice 
that  for  this  second-order  moment,  the  amount  of  noise 
in  the  simulation  data  has  increased  over  that  shown  in 
Figure  3.  This  is  due  to  the  fact  that  it  is  more  difficult 
to  obtain  higher-order  statistics  from  a  simulation,  often 
creating  more  ambiguity  in  the  data.  Despite  the  noise, 
however,  we  see  that  the  MEMC  mode!  agrees  well 
with  the  simulation. 

Figure  6  shows  the  mean  drop  axial  velocity 
profiles.  The  region  closest  to  the  injector  is  where  the 


Figure  6  Expected  drop  axial  velocity 

highest  velocity  slip  occurs.  There  are  a  number  of 
drops  passing  through  this  region  that  are  moving  faster 
than  the  mean  gas  velocity.  Due  to  their  lower  masses, 
most  of  those  drops  rapidly  decelerate,  reducing  the 
mean  velocity  towards  the  mean  gas  velocity.  At  the 
leading  edge,  however,  the  mean  velocity  increases.  To 
have  reached  the  leading  edge,  those  drops  had  to  be 
the  fastest  coming  out  of  he  injector,  and  because  they 
are  also  some  of  the  largest,  they  haven’t  been  affected 
by  drag  as  much  as  the  smaller  drops.  As  the  wave 
front  continues  to  move  forward,  those  large  drops  slow 
as  evidenced  by  the  decrease  in  the  mean  velocity  at  the 
front  as  time  goes  on. 

Figures  7  and  8  show  the  standard  deviations  of  the 
axial  and  transverse  velocity  distribution  functions.  The 
largest  standard  deviations  occur  near  the  injector 


Figure  7  Drop  axial  velocity  standard  deviation 


because  the  drops  have  not  been  influenced  by  the  gas 
phase  long  enough  to  have  approached  the  gas  velocity. 
Downstream,  the  standard  deviations  rapidly  decay 
until  we  reach  the  wave  front.  It  will  take  a  longer  time 
for  the  larger  drops  to  reach  the  gas  velocity,  but  as  the 
leading  edge  propagates,  the  standard  deviations 


Figure  8  Drop  transverse  velocity  standard  deviation 


become  smaller  and  smaller.  Comparisons  with  the 
Lagrangian  simulation  again  show  good  agreement. 
This  suggests  that  the  submodels  used  to  close  the 
cross-phase/same-component  velocity  moments  and 
term  appearing  in  the  turbulence  equations  representing 
the  effects  of  the  drops  do  a  satisfactory  job  in  allowing 
the  overall  model  to  predict  the  statistics  of  the  spray. 

Conclusions 

We  have  incorporated  the  effects  of  nonlinear 
momentum  exchange  between  the  phases  into  a  model 
which  describes  the  evolution  of  a  spray  flow  by 
solving  a  series  of  moment  equations  with  a  maximum 
entropy  model.  The  model  was  tested  for  a  quasi-one- 
dimensional  spray  flow.  Results  were  obtained  for 
several  quantities  of  interest,  including  the  expected 
drop  number  density,  expected  diameter,  expected  axial 
velocity,  and  various  second-order  moments.  For 
comparison,  a  Lagrangian  simulation  was  performed 
which  showed  good  agreement  with  the  MEMC  results. 
While  the  point  was  not  emphasized  in  the  present 
work,  we  determined  that  the  MEMC  model  requires 
less  computational  time  to  provide  similar  statistics  as  a 
simulation.  Finally,  we  showed  that  the  submodels 
used  to  describe  the  cross-phase  velocity  moments  and 
the  turbulence  modification  allow  the  overall  MEMC 
model  to  provide  reasonable  agreement  with  the 
simulation  data. 
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